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ABSTRACT 

Current explicit integration techniques in fluid dynamics are deeply limited by the 
Courant-Friedrichs-Lewy condition of the time step progression, based on the adopted 
spatial resolution coupled with the maximum value between the kinetic velocity or 
the signal transmission speed in the computational domain. Eulerian implicit inte- 
gration techniques, even though more time consuming, can allow to perform stable 
computational fluid dynamics paying the price of a relatively larger inaccuracy in 
the calculations, without suffering such a strict temporal limitation. In this paper, 
we present a simple and effective scheme to perform Free Lagrangian Smooth Particle 
Hydrodynamics (SPH) implicit integrations in Semi-Lagrangian approach without any 
Jacobian matrix inversion operations for viscous Navier-Stokes flows. Applications to 
SPH accretion disc simulation around a massive black hole (MBH) in a binary stellar 
system are shown, together with the comparison to the same results obtained accord- 
ing to the traditional explicit integration techniques. Some ID and 2D critical tests 
are also discussed to check the validity of the technique. 

Key words: accretion, accretion discs - conduction - convection - hydrodynamics: 
methods: numerical, N-body simulations - binaries: close - stars: novae, cataclysmic 
variables. 



1 INTRODUCTION 

A time step restriction is always necessary for time de- 
pendent calculations in computational fluid dynamics. Cur- 
rently, such restrictions are needed for mathematical sta- 
bility reasons in explicit calculations of partial differential 
equations (PDE), while they are necessary for accuracy con- 
siderations in implicit calculations. The integrated physical 
local property at time step level n + 1: A"^^ is a function 
of its previous values at time steps n, n — 1 etc., as well as 
of spatial derivatives of its spatial flux densities: dF{A)/dr, 
relative to the previous time steps for explicit calculations. 
Instead A"^^ is a function of these quantities also for the 
same n + 1 time step level for implicit calculations. 

For explicit calculations of the co mputational flow, the 
Cour ant-Friedrichs-Lewy condition (jCourant et al.l 1 19281 . 
1 19671 ) is imposed on those hyperbolic terms representing ad- 
vection in PDE (spatial derivatives of pressure or velocity), 
where the given Courant number C — VcAtcpL/^r ^ 1 
is generally of the order of 0.2 — 0.5, where Ar is the spa- 
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tial resolution, Vc is the maximum value between the local 
kinematic and the signal transmission velocities within the 
entire computational domain, and AtcpL is the Courant- 
Friedrichs-Lewy time step to be computed. 

Explicit integration techniques are widely adopted to 
solve equations of the fluid dynamics both in the Eule- 
rian formalism, where time and space derivatives refer both 
to local derivatives of the physical properties (d/dt and 
d/dr), according to the adopted spatial resolution length, 
and in the Lagrangian formalism, where the material or 
the convective deriva tive d/dt = d/dt + v ■ V characterizes 
the flow description (iFletcheij [l99ll : iHirschI Il997l : iLeVegud 
2002). Although the scientific literature is quite rich as for 
the implicit numerical integrations of PDE, implicit back- 
ward difference methods are also shown for ordinary dif- 
ferent ial equa tions (ODE) working with Lagrangian deriva- 
tives (jSewell 1988; Shampine 1994'). In spite of such a di- 
chotomy Ardoljian ct al. (1996) produced an implicit La- 
grangian method, based on a triangular mesh for calcula- 
tions of non stationary astrophysical processes, whose re- 
sults successfully compare with other Lagrangian explicit 
calculations, although some details in the flow do not com- 
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pare with those obtained via Eulerian formahsm. In Semi- 
Lagrangian techniques, the basic idea is to discretize the La- 
grangian derivative of the solution in time instead of the Eu- 
lerian derivative. Despite the full Lagrangian nature of the 
technique, a spatial grid is necessarily considered. Values of 
physical quantities of the system of differential equations are 
calculated through spatial interpolations at the grid points. 
Explicit-implicit integration techniques involve these spatial 
grid points and particle positions at previous time levels . 
Semi-Lagrangian schemes (|Robertlll969l : [Robert et aLlllQ?^ ) 
have been built up based on mixed explicit-implicit schemes, 
increasing the time step up to a factor of 6, paying a little 
additional cost in the time computing, with a limited degra- 
dation in the accuracy of solutions. However, the evaluation 
of the maximum stable time step still remains debated be- 
cause a time step lengthening of a factor of 6 appears much 
short er than seems necessary from considerations of accu- 
racy (|Robert|[l98ll ). These works of the 70's - 80's were 
developed in metereology for numerical weather forecasts, 
where the use of a longer time step is essential for efficiency. 
2D and 3D applications of integrations in Semi-Lagrangian 
approach have been discussed in Staniforth & Goto (1991) 
following a hierarchy of cases. 

The extension of Semi-Lagrangian method to the so- 
lution of viscous Navier-Stokes equations was done by 
IPirroneaul l|l982l ). showing the nonlinear stability of the in- 
tegration method, even as the viscosity decreases to zero, 
also s h owing an estimate of erro r s, lat er improved by ISiilil 
l|l988h : IXiu and Em KarniadakisI (|200ll) . 

Currently, SPH hydrodynamics works adopting an ex- 
plicit integration t echnique, being SP H a "Free Lagrangian 
particl e scheme" ( Whitehursd 1995[). Recent l y some au- 
thors llWhitehouse fc Bate! |2004| : ISusal l2006l : iRook et all 
I2OO7I : IPetkova fc Springelll2009l ) developed implicit integra- 
tion schemes working in SPH to solve some selected prob- 
lems dealing with the radiative transfer or with the heat 
transfer in the flow, adopting either the Crank-Nicolson, 
or the Runge-Kutta-Fehlberg numerical integrations, or the 
conjugate gradient method. The adoption of such meth- 
ods involves the handling of some time-expensiv e and mem- 
ory c onsuming Jacobian matrices. Alternatively l|Viau et al.l 
I2006D . some iterative variational techniques are also used 
to find zeros for the analytical equation coming out from 
the energy equation written in finite terms and including 
the diffusive terms, within an assigned tolerance error, set- 
ting the "left" and "right" boundary values for the ther- 
mal energy per unit mass e, finding the median value us- 
ing the Van Wii ingarden-Dekker-Brent bisection method 
l|Press et al.lll992l ). However, the introduction of such mixed 
procedures involves deviations from the original pure parti- 
cle hydrodynamics. 

In this paper we present a simple implicit mathematical 
technique able to perform Semi-Lagrangian explicit-implicit 
integrations both of the Euler and of the Navier-Stokes equa- 
tions of the fluid dynamics, respecting the SPH, or the SPH- 
derived schemes, without any Jacobian matrix handling. 

The implicit integration procedure in our scheme is nat- 
urally based on an iterative integration scheme. However, we 
pay attention to the solution of the entire system of equa- 
tions for all SPH particles at the same time instant without 
using interpol ation techniques like other Sem i-Lagrangian 
schemes need (|Robertlll969l : lRobert et al.lll972l ). limiting the 



procedure to the i mplicit solution of a sing le advection equa- 
tion even in 3D (jStaniforth fc C6tg|l99ll ). Instead, we use 
SPH interpolations to perform implicit iterative integrations 
on the entire system of equations for all SPH particles within 
the computational domain. 

An implicit integration scheme promises to be more ef- 
ficient in some situations, for example when large bulk ve- 
locities are present, or when a system is close to hydrostatic 
equilibrium, or for very low Mach number flows, for string 
diffusion, or for stiff source terms, etc.. In such cases the 
Courant-Friedrichs-Lewy time step constraint that needs to 
be obeyed in explicit time integration schemes can be very 
restrictive, and may yield an explicit time step much shorter 
than in principle needed to follow the real dynamics of the 
system. 

The next sections of this work show how a Semi- 
Lagrangian explicit-implicit technique can be applied to 
SPH. To do this, we brieffy recall how SPH works in §2. 
Then, we show how implicit integrations schemes work, how 
the Lagrangian Euler or Navier-Stokes equations have to be 
rewritten, how they are translated in the flnite-difference 
schemes and how a suitable implicit time step advance- 
ment well beyond the Courant-Friedrichs-Lewy limit is cho- 
sen (§3). Successful applications are here also reported for 
the case of ID and 2D blast waves SPH modelling, as well 
as for the 2D gas shockless radial viscous transport (§4). Fi- 
nally, astrophysical applications are here also presented (§5) 
on the study of a 3D an accretion disc around a MBH in 
a microquasar, both adopting the implicit and the explicit 
integration procedures. Some advantages and some disad- 
vantages of the implicit techniques, compared with the ex- 
plicit schemes are also discussed. Both explicit and implicit 
accuracies, applied to SPH, are also discussed. 

In a physically viscous accretion process onto a MBH in 
a microquasar, the viscous conversion of mechanical energy 
into heat decreases the Mach number to values so low that 
the CFL explicit time step could be too short for practical 
purposes. 



2 FLUID DYNAMICS EQUATIONS AND 
THEIR SPH FORMULATION 

The relevant equations to our model or viscous gas hydro- 
dynamics are: 



^ + pV • u = 



continuity equation(l) 



at p p 

Navier-Stokes momentum equation 



energy equation(2) 
equation of state(3) 

kinematic equation. (4) 



dr 



dt 
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d/dt stands for the Lagrangian derivative, p is the gas 
density, e is the thermal energy per unit mass, / is an ex- 
ternal force field, p is the ideal gas pressure, here generally 
expressed as a function of local properties, determined by its 
equation of state (EoS). The adiabatic index 7 has the mean- 
ing of a numerical parameter whose value lies in the range 
between 1 and 5/3, in principle, r is the viscous stress ten- 
sor, whose presence modifies the Euler equations for a non 
viscous fiuid dynamics in the viscous Navier- Stokes equa- 
tions. Notice that the inclusion of the field terms (inertial 
or non inertial) does not affect the mathematical scheme re- 
garding the theme of this paper. Sometimes these terms are 
absent and sometimes they are present, according the sim- 
ulation we are considering. Here they are written only for 
universality reasons. 

The SP H method is a Free Lagrangian scheme 
l|Whitehurstl [1995 ) that discretizes the fluid into moving 
interacting and i nterp o lating domains called " part icles" 
ljMonaghanlll985l . Il992l : iMonaghan fc Lattanziol Il985l l. AU 
particles move according to pressure and body forces. The 
method makes use of a Kernel W useful to smoothing inter- 
polate a physical quantity A{r) related to a gas particle at 
position r according to: 

'A{r)= / A{r')W{r,r' ,h)dr' . (5) 

J D 

W{r,r' ,h), the interpolation Kernel, is a continuous 
function - or two connecting continuously differentiable func- 
tions even at the connecting point - defined in the spatial 
range 2h, whose limit for /i — > is the Dirac delta distribu- 
tion function. All physical quantities are described as exten- 
sive properties smoothly distributed in space and computed 
by interpolation at r. In SPH terms we write: 



^A. 



^ — ' n, ^ — ' n,- 



(6) 



where the sum is extended to all particles included 
within the domain D, rij = pj/rrij is the number density 
relative to the jth particle. W{ri,rj,h) is the adopted in- 
terpolation Kernel whose value is determined by the relative 
distance between particles i and j. J W{ri,rj,h)d^r' — 1, 
that is: W{ri,rj,h)/nj = 1. 

In SPH conversion of mathematical equations there are 
two principles embedded. Each SPH particle is an extended, 
spherically symmetric domain where any physical quantity 
/ has a density profile fW{ri,rj,h) = fW{\ri — rj\^K) = 
/W^d^ijl, /i). Besides, the fiuid quantity / at the position of 
each SPH particle could be interpreted by filtering the par- 
ticle data for /(r) with a single windowing function whose 
width is h. So doing, fiuid data are considered isotropically 
smoothed all around each particle along a length scale h. 
Therefore, according to these two concepts, the SPH value 
of the physical quantity / is both the overlapping of ex- 
tended profiles of all particles and the overlapping of the 
closest smooth density profiles of /. This means that the 
compactness of the Kernel shape gives the principal con- 
tribution to the interpolation summation to each particle 
by itself and by its closest neighbours. In both approaches 
the mass is globally conserved in so far as the total particle 
number is constant. 

In SPH formalism, equations (2) and (3) take the form: 



dvi 
~dt 



N 



Pi 



N 



Pi 



+ 



(7) 



dl^' 



N 



mi 



ViWij + Qi ■ Vi + 



N 

E 



mj \ rivi — ^2 h rivj^^-^^ ) ■ ^iWij (8) 



(Ti ■ Vi 



where rrij is the mass of jth particle and p* = Pi+ dis- 
sipation pressure term. Ei = (e^ -f ^v^)- The viscous stress 
tensor Tap includes the positive first and second viscosity 
coefficients r]v and C,v which are velocity independent and 
describe shear and tangential viscosity stresses ijiv), and 
compressibility stresses 



Tap = Vv'^aP + CvV ■ V 

where the shear 
dva , dvp 2 



(9) 



(10) 



In these equations a and P are spatial indexes while 
tensors are written in bold characters. For the sake of sim- 
plicity we assume C^v ~ 0, however our code allows us also 
different choices. Defining 



,^ \ ^ rrijVjia dWij 

ViaP — > 5 



(11) 



as the SPH formulation of dva/dxp, the SPH equivalent 
of the shear is: 



CiaP — ViaP + ViBa — —5apVi^-y 



(12) 



A full justi fication of this SPH formalism can be found 
in lFlebbe et~ ( 1994a, b). 

In this scheme the continuity equation takes the form: 



JV 

dpi _ \ ^ 

~dt 



/ rrijVij ■ ViWij 



or, as we adopt, it can be written as: 



JV 

Pi = ^mjWij 



(13) 



(14) 



which identifies the natural space interpolation of par- 
ticle densities according to equation (7). 

In such a conversion, the physical mass density of the 
ith SPH particle is either physically represented as pi = 
miTii, or numerically represented as pi = pj /rijWij = 
. rrijWij. In the same fashion, the particle numerical den- 
sity is either Ui = pi /rui , or = Wij . 

A necessary convergence is needed, because the two ex- 
pressions for Pi or for Ui coincide only in the case of equal 
mass SPH particles, when rrii — rrij. 
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In addition, looking at the SPH particle masses, either 
rrii = pi/ui, or m; = J^j ^i/'^i^ii = Ylj Pi/'^'j^iJ' which 
equal with each other only when m = nj . 

In some circumstances SPH particles could have differ- 
ent physical masses. In this case, the concept of SPH particle 
mass could be without any physical meaning whenever local 
densities are different: pi 7^ pj. 

To resolve this incongruity, we modify the SPH formu- 
lation (15) according to one of these choices: 

either 



m — = — . mjWii 
pi = miUi = Ylij "^i^ii-, 

as we did in this paper, or 



(15) 



pi = m,n, = Hi J2j -^Wij 
or 

Pi — miUi — rniJ2i 



(16) 



(17) 



The first of these p ~ n ~ m correlations (eqs. 16) is 
written considering the SPH interpolation integral (eqs. 6, 
7) applied only to the mass density p, not to m or to n. The 
second correlation (eqs. 17) is obtained considering eqs. (6, 
7) applied only to the mass m, not to p or to n; instead the 
third form (eqs. 18) comes out applying eqs. (6, 7) only to 
the numerical density n, not to p or to m. 

Being p, n and m physically correlated as p = mn, 
the SPH interpolation integral (eqs. 6, 7) cannot be used 
at the same time for all these three quantities, otherwise 
differences among smoothed to unsmoothed quantities could 
be relevant. 

Instead, according to one of eqs. 16, 17 or 18 choices, the 
full convergence between the SPH transformations and the 
full physical meaning of mass and densities are respected. 

In its original SPH formulation, standing the ideal equa- 
tion of state (EoS) in the form: 



(7 - l)pe 



perfect gas equation, (18) 



dissipation in p* is given by an artificial viscosity term. 
This dissipation, together with an appropriate numerical 
thermal d iffusion con t ribut i on oc ( Uj — Uj)csii / Pa, where 
Ui = piti iMonaghanI (|l985l . ll992D : iMonagh an fc Lattanzid 
l|l985h . included in dt/dt, reduce shock fluctuations. The ar- 
tificial viscosity term is given by: 



Vij = asPHP-ij + PsPHPij, 



where 



^2— 



P-ij — 







if Vi 



< 



(19) 



(20) 



otherwise 



being Csi the sound speed of the ith particle, Vij — 
n - Vj, Vij = Vi - Vj, < h^, asPH ~ 1, PsPH ~ 2, 
Csij = 0.5(csi -I- Csj) and pij = 0.5(pi + pj). These asPH 
and PspH parameters of the order of the unity are usually 



adopted to damp oscillations past high Mach number shock 
front s developed by non-linear instabilities (|Boris fc Bookl 
Il97i ^. These asPH an d Psph values were also adopted by 
iLattanzio et al.l (Il985l). Smaller a sPH and Psph values, as 
adopted by Meglicki et al.l (|l993l ). for developing more tur- 
bulence. In the physically inviscid SPH gas dynamics, an- 
gular momentum transport is mainly due to the artificial 
viscosity included in the pressure terms as: 



2i 

Pj \P^ Pj 



(21) 



where p is the intrinsic gas pressure. 

However, looking at a physical ori gin of the numeric al 
dissipation term, inclu ded in the EoS (|Lanzafamel I2010al l3: 
iLanzafame et al]|201ll ). gas pressure can be expressed as: 



C- 



3cs 



(«7:) 



(22) 



(23) 



P 2 

P = -Cs 

1 

where 
^ 1 

C = — arccot 

where i? 3> 1. -R is a large number describing how much 
the flow description corresponds to that of an ideal gas: 
R ~ A/d, being A oc p~^^^ the molecular mean free path, 
and being d the mean linear dimension of gas molecules, vp 
is the impact relative velocity component. The physical dis- 
sipation, expressed by the two further terms in eq. (23) (the 
linear and the quadratic terms in V • d) of the reformulated 
EoS, better treats both shocks and shear flows, even in a 
Lagrangian description. Their inclusion substitutes artificial 
viscosity terms and does not represent a physical turbulent 
viscous contribution, but the physical dissipation coming out 
because eq. (19) of the EoS should strictl y be applied only to 
quasi -static processes (|Lanzafamell201o"al 3: ILanzafame et al.l 

[2oi3)- 



3 CONCEPTS ON EXPLICIT AND IMPLICIT 
INTEGRATION SCHEMES: APPLICATION 
TO SPH 

3.1 Formulations on explicit and implicit 

generalized three-level integration schemes 

Given a physical property A, the mathematical conversion 
of its time derivative dA/dt and of its first and second 
space derivatives: dA/dr and A/dr^ is: dA/dt — i- {A 



n + l 
k 



Al)/l^t, dA/dr 0.5(A 
0.5(A 



A^_i)/Arfe and d^A/dr^ 



fc+i 



- 2yl^ + Afc_i)/Arfc, respectively for explicit tech- 



niques, while it is dA/dt -5> (^4^+^ 
Q.^{Al+l~Alt\)/^rk and d'^A/dr'^ 
Al" 



A^)/At, dA/dr 
■ 0.5(A- " 



)/Arfc, respectively for implicit techniques (|Fletcheil 



I1991I : lHirsc"h| [r997: LeVeque _2002). In such expressions, 
represents the temporal level, while k represents the spatial 
grid index, according to the versus of the reference frame. 
At = - and Arfe = ru+i - rt = rk - r^-i = 

0.5(rfc+i — rfc_i). Of course other higher order sche mes exist, 
where more time and space leve ls are considered ( |Fletcherl 
I1991I : lHirschlll997l : iLeVeaudlioO^ ) , which we do not consider 

for the sake of simplicity. 

As a useful generalization (|Fletcher|[T99ll '). the hypo- 
thetical equation: 
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dt dr dr^ 



c = 



(24) 



can be converted as: 



+ 

a (1 - 



- A- 



+ K- 



+ 



= 0. 



(25) 



+ 



Whenever the pair A = and k = 0, such a general 
expression yields the typical explicit two points forward cen- 
tred integration scheme (2FCS) for PDE. Instead, for A = 
and At = 1 we have a simple full implicit centred three points 
integration scheme (3FICS). Furthermore, for A = 0.5 and 
K = 1 we have a linearized full implicit three points tech- 
nique (3LFI), while for A = and k = 0.5 we obtain the 
well known Crank-Nicholson implicit scheme (|FletcheJll99ll : 
lHirschlll997l : lLeVeaue|[20o3 ). 

The truncation error jFletchej [T99ll ') is a function oc 
Arf{d^A/dr^,d^A/dr^) for upwind 2FCS, as well as for 
3FICS schemes, while it is oc Ar'^ f(d^A/dr^ ,d^A/dr'^) for 
the other last two implicit techniques. 

Although the expressions (25) and (26) refer to the nu- 
merical value of a non-smoothed or intrinsic physical prop- 
erty A, we assume that it also holds for its smoothed SPH 
evaluation. Hence we can rewrite the same expressions us- 
ing indifferently A instead of A in eqs. (6) and (7). Since 
we integrate both explicitly and implicitly only smoothed 
physical quantities, throughout the rest of the paper we will 
use mathematical symbols without any superscript line for 
a simpler reading. 



(k-ln + l) (kji+1) (k+l,n + l) 



(l,f) 




(l.f-) 



Figure 1. Schematic plot showing the position of the ith SPH 
particle at time t"-^,t" and At time the X,Y,Z spatial 

grid is shown, where physical properties are interpolated, used to 
perform the implicit 3FICS scheme, together with the information 
relative to the same particle at previous times. A stencil diagram 
referring to the ID implicit integration scheme is also reported on 
the left of the same picture showing the spatial k and temporal n 
indexes. In 3D, a distorted stencil scheme is used in like manner 
as in ID. 



3.2 SPH in a Semi-Lagrangian explicit-implicit 
generalized three-level integration scheme 

Hydrodynamics in the nonlinear Free Lagrangian SPH ap- 
proach is currently performed in predictor-corrector explicit 
schemes, starting from some initial values at time t = 0. In 
the "Leapfrog" scheme, the equations for space and velocity 
advancement can be written as: 

Tk = r^+v^ 'At (26) 

n+1/2 n— 1/2 . n a , /otX 

«r-.fe = "r,fc +ar,k^t (27) 

that can be manipulated into a form which writes par- 
ticle velocity at integer steps as 

n+l n , n \ , , ^ n \ ,2 /(-\n\ 

rk = r-fc + fr,feAt -I- -ar,fcAt (28) 

Wr,fc = Wr.fc + gK.fc +"r,fc )At. (29) 

In this second expression, since particle acceleration 
a depends on v, it is required an implicit integration for 
the second equation. In the case of a "Leapfrog" scheme, 
an "evaluator" phase in the computational scheme needs 
to be interposed between the two integration procedures, 
where time derivatives of the various physical quantities are 
computed. For this reason, this scheme is a so called PEC 



method, where a Predictor-Evaluator-Corrector procedure 
is followed by the updating of all integrated values. 

Iterative Runge-Kutta methods are also used, both ex- 
plicit as well as implicit. In such schemes, the integrated 
value of the physical property A^'^^ — A^ + SAt, where S is 
an estimated weighted average of slopes from the beginning, 
through some midpoints, toward the end of the time inter- 
val. Despite more general than explicit methods, implicit 
Runge-Kutta schemes are more complicated and dependent 
on the specific problem. Those Runge-Kutta methods that 
are diagonally implicit, show a strong stability allowing a 
significant increase in the time step li mit, compared with 
the explicit methods of the same order (jVisbal fc Gaitondd 
1 19991 : 1 Ketcheson et aLll2009l ). Here, we do not discuss in de- 
tail this specific complex argument, where often some Jaco- 
bian matrices need to be inverted. However, even the simple 
backward Euler (A = 0, k = 1), or the Crank-Nicholson 
(A = 0, K = 0.5) methods belong to this category. 

The explicit multistep Adams-Bashforth-Moulton 
PECE explicit integration scheme can also be adopted in 
SPH, where either the Adams-Bashforth 
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A2+^ = A?.+ At- 



at 



= + ^ 3 



1 , At / oS^fc 



at 



Al+^ + # I 23:^ - 16^9^ + 5 



at 



£££ 

at 



"^fc +24!°° at °^ at 



+ 



(30) 



eA^_ aA£ 
at ^ at 



■ 2774 



aA" 



-+ 



2616 



at 



1274 



aA 



— +251^ , 



or the Adams-Moulton 



aA" 



at 



+ 



J^+Z _ _j_ At I ^ 



Al+''^ = ^fc^^ + # I 9 



at """at 



SA" 



at 



+ 8 



8A" 



££2 

at 



aA" 



+ 19 



SA" 



(31) 



n+3 



+ 



251- 



at 



646^- 264^ + 106^-19 



at 



are used. To apply either the Adams-Bashforth or the 
Adams-Moulton techniques, up to the wished precision, pre- 
vious derivatives for the same flow elements need to be con- 
served. Besides, a further evaluator procedure is considered 
at the end of the predictor-corrector integration scheme, to 
be a PECE (not a PEC) technique. 

What has been discussed is necessary to under- 
stand how it is possible to build up a Semi-Lagrangian 
explicit-impHcit technique for SPH. In our Semi-Lagrangian 
explicit-implicit SPH integration approach, the first ex- 
plicit integration scheme is necessarily the same that 
we currently use (Leapfrog, Adams-Bashforth-Moulton 
schemes, etc.), without strictly taking into account of the 
Courant-Friedrichs-Lewy limit for the time step. The adop- 
tion of an explicit-implicit backward and forward proce- 
dure has also been used by some ot h ers l| Miranda et ahl 
19891: iGravouil fc Comberscurd I2001I : iMaiid et al.l I2OO6I : 



Ying et al.ll2008l : llsmail et al.l bOOgfT although some recent 
pure implicit sche mes exist (S ahin fc Qwii3 |2003| . l2006l : 
iMosgueda fc Ahma dizadch 2010), under some strict condi- 
tions. We use indifferently either a Leapfrog or an Adams- 
Bashforth-Moulton scheme as explicit procedure to calcu- 
late integrated values to the time level n + 1 for all SPH 
particles. Then, we can use iteratively the expression (26) 



to repeat, only implicitly, integrations adopting A = 0.5 
and K, = 1, without any mathematical operation of Ja- 
cobian matrix inversion. The total number of implicit it- 
erations is of course arbitrary. The first explicit integra- 
tion step in the Semi-Lagrangian explicit-implicit technique 
works only once at the 1st iteration. From the 2nd iter- 
ation onward, the entire integration process is fully im- 
plicit. However, 2 or 3 implicit iterations at most are usually 
enough to get a convergence (even better in SPH) because 
the variation between the implicit solution and the previ- 
ous solution in the iterative integration procedure quickly 
reduces toward very small relative differences Au/u (being 
u = p,t,vx,VY,vz, avoiding any divergence in the solu- 
tion (|Sulilll988l : [Xiu and Em Karniadakia,200L '). Notice that 
these relative differences refer to relative variations of the 
generic integrated physical quantity u by the progression of 
iterative calculations of computed integrals, not to relative 
errors in the implicit integrations. Relative errors refer to 
relative discrepancies from the exact analytical solution, if 
it is known, or to general aspects referred to the entire flow 
through conservation laws. A number of implicit iterations 
larger than 3 is of course possible, but practicable in so far 
as the implicit iteration technique is competitive against the 
current explicit integration methods. On the contrary, not 
only results would be less accurate, but also the entire pro- 
cedure would be more time consumi n g. In t his sense, the 
Semi- Lagrangian schemes of iRobertI (|l969l ): [Robert et al.l 
(|1972|) . whose method works for a time lengthening up to 
a factor of 6, should be not convenient if implicit iterations 
are more than 4 — 5. A quick convergence of numerical re- 
sults in SPH is ensured by the fact that dissipation (even 
artificial) always works since the first explicit time step of 
the iterative cycle so that solutions normally relax without 
any oscillation or with very small damped oscillations. A 
larger number of iterations also produces a cumulative ef- 
fect both on dissipation and on numerical errors, cycle by 
cycle, which is normally a not wished result. A control on the 
solution accuracy is normally possible in so far as only one 
SPH equation is implicitly integrated (as some authors do 
only for the energy equation). But, this is practically impos- 
sible when solving implicitly the entire system of equations, 
because how to control at the same time so many Semi- 
Lagrangian solutions as many are the scalar equations for 
a multitude of grid points is not treated in the mathemati- 
cal literature, without any warranty for an improvement on 
other solutions. As an example, it could be possible that 
the implicit solution for the density for one SPH particle 
satisfies the arbitrary criterion that its relative variation is 
limited within the assigned tolerance as Ap/p < 10~^ after 
3 implicit iterative cycles, but at the same time that the 
other relative variations for Ae/e, Avx/vx, etc. still do not 
because they are greater than the established tolerance. The 
iterative cycle could either terminate, accepting the relative 
variations as they are at the moment, or it could continue. 
But each further cycle will destroy the consistency of the 
solution for the density which relaxes toward its numerical 
value within 3 cycles of implicit iterations. From the 4th 
implicit cycle ahead the solution for p could loose accuracy 
more and more. Even though normally a local instability af- 
fecting one solution (explicit or implicit) also affects other 
integrated solutions, situations where discrepancies occur in 
the relative variations cannot be excluded. These difficulties 
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normally increase as longer is At; compared with AtcFL 
because of both error and dissipation accumulation. When- 
ever Ati/ AtcFL < 6, the degradati on in the accura cy is 
normally negligible (|Robert. 1969; Robert et alJIlOTd ). but 
for larger time ratios it is not. This means that, especially 
for 10 < At/ AtcFL < 20, the number of iterations in the 
implicit procedure cannot be too low because the solution 
could be affected by the instabilities coming from the first 
explicit integration, but also it means at the same time that 
the number of iterations cannot be too high because of the 
long time lost and because of the inaccuracy corruption on 
the solution. A relaxation in the SPIf solution normally hap- 
pens. The density of energy U = ^pv^ + pe + pf ■ r or the 
density of enthalpy H = U+p include all physical quantities 
p, V, e, However, the propagation of relative discrepancies 
Au/u for both p, e, v^, Vy and Vz cannot be managed cycle 
by cycle. What is here discussed clearly shows that the ar- 
gument of solution consistency for a system of equations is 
quite complex for only one particle. However each iterative 
cycle of implicit integrations needs to involve a multitude 
of SPH particles at the same time. Relative variation for 
energy or enthalpy densities for the totality of particles as 

^Ui/ or as AH^/ < an arbitrary tol- 

erance (e.g. 10~^) is a valid general criterion for the evalua- 
tion of the consistency of the implicit integrated solutions at 
the same time level n + 1 between two consecutive integra- 
tion cycles. These global conditions do not ensure the con- 
sistency of local implicit solutions for the single particle even 
for explicit integrations within the CFL condition. This, in 
particular happens whenever free edges of the flow are con- 
sidered, without any Dirichlet or Von Neuman boundary 
conditions. Indeed, two different SPH particles could have 
very different values for AUi/Ui or for AHi/Hi, especially 
at the free edge of the computational domain. Hence, even if 
this global criterion is not satisfied at the 4th iterative cycle 
(1 explicit plus 3 implicit), we decide to interrupt implicit it- 
erations, proceeding to the next time step. Thus we privilege 
the criterion that the integration for each implicit iteration 
cannot exceed the 4th cycle (the third implicit cycle), oth- 
erwise the whole implicit scheme would be not competitive 
compared with the explicit integrations, because too time 
consuming. Since this paper is a study on the application 
of a Semi-Lagrangian integration technique to SPH, even 
the same comparison with SPH explicit integration solution 
should answer to the question whether the Semi-Lagrangian 
explicit-implicit integration consistency of solutions is ac- 
ceptable or not. 

Our Semi-Lagrangian explicit-implicit integration tech- 
nique on SPH costs about 10% more time for each iterative 
implicit cycle (from the 2nd cycle onward) compared with a 
simple Leapfrog explicit scheme as used in SPH. This is due 
to the fact that the most expensive calculation is that rela- 
tive only to the SPH interpolations of the two ends of each 
spatial stencil at time level n -I- 1 on the ith particle, fortu- 
nately restricted only to the same known local neighbours. 
This means that it is not necessary further computational 
time in searching for neighbours to the ith particle. At the 
same time, it is also necessary a limitation of the implicit 
dissipation accumulation. Truncation errors consist, typi- 
cally, of succes s ive higher even and odd ordered derivative s 
l|Fletchedll99ll : lLeVeaue|[l99^ : lHirschlll997l : lLeVeaue|[200^ ) 
acting as a dissipation on the integrated solution. Therefore, 



each implicit iteration yields a solution including a dissipa- 
tion, involving further stability, but also a larger inaccuracy, 
as it is normal for implicit integration calculations. It is nor- 
mal that the implicit integrations are anyway less accurate 
than the explicit ones. However a much shorter duration 
of calculations should well compensate this disadvantage if 
consistency of solutions is satisfactory. 

Notice that in expression (26) while the r subscript in- 
dicates the component of a vector determined by the orthog- 
onal projection, the k subscript is an index of position along 
a ID spatial mesh at time (n + 1), where SPH interpola- 
tion are needed. We adopt this scheme in 3D along three 
arbitrary distinct orthogonal axes (e.g. axes parallel to the 
X, Y, Z directions), assuming that the real SPH particle 
stays at the centre of the spatial-temporal (k, n+l) ID line 
parallel to the chosen arbitrary axis, time by time (Fig. 1), 
whose spatial grid index k runs along it. The grid index k 
identifies the ith particle along the ID computational arbi- 
trary spatial line (here X, Y, or Z) in 3LFI. Indexes fc -I- 1 
and fc — 1 refer to the interpolation spatial points "ahead" 
and "behind" the same ith SPH particle along such an axis. 
The necessity of as many orthogonal axes as many are the 
dimensions of the flow is necessary to solve discontinuities in 
the flow whenever they are not simple linear fronts. So do- 
ing, we assume that both the Eulerian and the Lagrangian 
spatial derivatives are expressed in the same fashion in the 
finite differences scheme and that the correspondence be- 
tween the temporal Lagrangian derivative d/dt and its cor- 
responding terms in the first row of eq. (26) still holds as for 
the Eulerian schemes. This can be done because such finite 
SPH discretizations refer to the moving Lagrangian particle, 
not to a spatial grid cell. On the left margin of Fig. 1, it is 
shown the space-time regular grid th at is typically consid- 
ered for ID finite di fference schemes feahin fc Owensll2003l . 
l2006l: iMosgueda fc A hmadizadch 2010) for a direct compar- 
ison to the distorted stencil scheme shown in 3D in the same 
picture. In reality, although the (fc, n) space-time mesh term 
is technically correctly used, we do not use a real space-time 
mesh as for a pure Eulerian scheme in the 3D space, but sim- 
ply some stencils connecting the (fc -I- 1, n + 1), (fc, n + 1), 
(fc — 1, n + 1) positions at time n+1 with the earlier (fc, n), 
(fc, n — 1) positions of the ith SPH particle (Fig. 1). 

In the SPH approach, since smoothed solutions A are 
implicitly integrated instead of A (§2, §3.1), as an indi- 
rect benefit, the convergence of solutions toward their val- 
ues should be effective since the second implicit iteration. 
This normally does not suddenly happen for non smoothed 
physical solutions. Therefore, the Semi-Lagrangian explicit- 
implicit integration approach well combines with SPH in so 
far as the two pseudo-particles ahead (fc -I- 1, n -I- 1) and be- 
hind (fc — l,n + 1) are located within 2h from the central 
real particle position at (fc, n+ 1). This is the reason why we 
consider the spatial lengths from (fc — 1, n -I- 1) to (fc, n+1) 
and from (fc, n -I- 1) to (fc -I- 1, n -|- 1) positions equal to h. 

The characteristic of Semi-Lagrangian explicit-implicit 
techniques is that the full Lagrangian description is con- 
served. The first step of the integration process is always 
explicit (as it is always necessary as for the full implicit 
techniques). Only physical values at the previous two time 
steps for each moving Lagrangian particle are needed. But, 
and this is the most important thing, only the physical infor- 
mation calculated at two adjacent fictional equally spaced 
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pseudo particles along a stencil, where the real particle stays 
at its centre, in the 3D space are necessary. Hence a few 
space-time mesh points are involved for each iteration on 
each particle. This has the advantage that in the Semi- 
Lagrangian explicit-implicit integration process, heavy ma- 
trix inversion operations are avoided. Pure implicit integra- 
tions, using 3D grids axe indeed memory and time consum- 
ing if a 3D simulation is considered, whose spatial resolution 
length is very short. The results are comparable, in princi- 
ple, but the Scmi-Lagrangian approach is technically faster 
and it requires a smaller computer memory. 

In the conversion of the physical Euler or Navier-Stokes 
system of equations in its discretized form using eq. (26), it 
is essential that dA/dt oc dF{A) / dr, where F(A) is the flux 
density of A. If dA/dt oc Adv/dr for example, some algebraic 
operations are needed to put A inside the spatial derivative, 
within the routines regarding only the implicit integration 
scheme. 

According to this strategy, adopting (A = 0.5, k = 1) 
the continuity equation, written as: 







(32) 



in finite difference terms is converted as: 



At 



2 At 
1 



2Arfc 



[Pk+l^r. 



2Ark 

n+l 
■ ■ 



{plXl 



Pfc^i<l-i)=0, (33) 



where r is the Cartesian component of the velocity vec- 
tor. Taking into account that A; is a ID grid number along 
each X, Y, Z arbitrary directions time by time, we have 3x3 
equations to be solved to get as many p^'^^ values to be av- 
eraged. In this last numerical expression the ID grid number 
k, shown as a subscript, indicates the two ID positions fc — 1 
and fc -|- 1 where the physical property is SPH interpolated, 
being fc the index pointing the real particle position, along 
the line connecting fc — 1 to fc + 1. In this Somi-Lagrangian 
explicit-implicit technique, this interpretation is applied also 
to any other scalar quantity or scalar vectorial component. 
Also the conversion of the divergence term in this numerical 
expression must be interpreted in like manner as an ID con- 
version. This is the reason because we get 3x3 equations to 
be solved. Even though we use the three X, Y, Z directions 
as preferred lines to work as fc directions, the choice of other 
arbitrary orthogonal directions does not affect the final re- 
sult in so far as the adopted SPH Kernel is isotropic as it is 
currently used in SPH. 

In some particular cases, it is better to rewrite the same 
continuity equation as: 



dlnp 
dt 



+ V ■v = 0. 



(34) 



In this case, its finite differences conversion is: 



3lnp^+i-lnp^ llnp^-lnpr 



At 



At 



- + 



^ fj.i+i _ ,,"+1 ^ - 



(35) 



As far as the energy equation is concerned, written as: 
■[(p + pe)v-VT\-vV{p + pe) = (36) 



we have: 



q n+l n+l _ .n,n -, „n n 
^ Pk Pk^k J- Pk^k 



„n — 1 n—1 



At 



At 



2^ [(P + P^)fe + l "rl+l -iP+ P^)k-1 «r,l-lj 
2Ark ■' V r.s.fc-Hl '•.fc-Hl ^r,s,k-l"r,k-l) 



2Ark 



[{p + perktl-iP + P^rk-l]=0, 



(37) 



where s = X,Y, Z. 

As far as the momentum equation is concerned, we 
postpone the inclusion of all explicit external field contribu- 
tions (e.g. gravitational, electric, magnetic etc.) after having 
computed the velocity due to the thermodynamics alone. 
Such contributions can be easily added to the integration as 
AtV^g^^^ f.. According to the same concepts, the momen- 
tum equation, without any external field contribution 



dVr _|_ 1 

dt p 



dp 



dr 

is converted as 



((V-T).r) 



= 



(38) 



3 <fe ' - <fc 1 'I'r.k - '-".k ' , Pk^t - Pk-l 



At 



At 



+ 



2pl+^Ar 



2pl+^Ark ^ ^'^"'^''+' " - 0- (39) 

Instead, if the same momentum equation is written as: 
^ - V ■ V(pw) -I- V • {pvv - r) Vp = 0, (40) 
its conversion in finite difference terms is: 



SPk - PkVr,k lPkVr,k-Pk Vr,k 



At 



At 



Ck' 



n + l _ n-H 

\Pk + l"r,k + l Pk-l"r,k-l) ^ 9> ^ 



2Ark v^^+i-'-.'^+i ,-^-,~r,^-,j . 

/ n+l n + l n + l _ n + l n + l n + l \ 
\Pk+l''^r,k+l''^s,k+l Pk-l'^r,k-l'^s,k-ll 



2Aru ^ 



^5^«;i+i-<;U)=0. (41) 



2Ark 



At last, as far as the Lagrangian position updating equa- 
tion dr/dt = V, it could be solved for each r = X,Y,Z 
component either explicitly as: 



~k 



: rl + Vr,kAt 



or implicitly as: 

..n+l ^n 1 ^n — 1 

2 ~At 



3 T'fc^'" ~ T'k 1 ''fc ~ r'k ..n+l 



At 



'^k,r 



: 0. 



(42) 



(43) 



Throughout these algebraic expressions, some kinetic 
and thermodynamic quantities have to be computed at (fc + 
1), (n-|-l) and at (fc — 1), (n-|-l) points in the spacetime grid. 
Adding the i subscript, referring to the ith SPH particle, we 
get: 
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and 



(44) 



(45) 



where As"^^ is a vector along the direction joining the 
space-time interpolating grid points (fc — 1), (n + 1) to {k + 
1), (n+l). Thus, the SPH technique is again used to compute 
the spatial gradients V Ai (e.g. gradients of density, pressure 
including its dissipation terms and so on) to the ith SPH 
particle within the scalar product in the 2nd term in eqs. 
(45, 46). For the sake of simplicity, even though |As"^^| is 



hi, whatever is the 



an arbitrary length, we adopt |As"^^| 
time level n. 

Once p^+\ e^+\ (r = X,Y,Z) are given at 

time n+l for the ith SPH particle for a = 1,2,3 in 



3D, the final values p" 



are calculated as: 



^n+i ^ y^^^y A"+79 where , where the two 

I ' 'X,Y,7 ' ' ^ / i,s ' 

summations take into account of the three flow components 
for each of the three directions (3x3). 

In this notation, the s index refers to the direction of 
the ID arbitrary computational spatial line, while the r sub- 
script refers to vectorial components. In our notation, s fol- 
lows r, but it could be different, in principle. Therefore, as 
an example, we compute p^+\ pl+\ p-+\ <+\ 6^,+^ 

that means the density and the thermal energy per unit mass 
on the ith particle, at time n-f 1, along the x, y, z lines (along 
the 1,2,3 arbitrary lines). Moreover, v^W, vl\\, v"+l, 

^T.' <,ti' <r.' <ti the r = x,y,z 

velocity components on the ith particle, at time n+l, along 
the s = x,y,z lines (along the 1, 2, 3 arbitrary lines). 

The strategy adopted in this explicit-implicit Semi- 
Lagrangian scheme is similar, for some aspects to that 
adopted to solve ordinary differential equations (ODE), 
even up to the second order in some e xplicit and implicit 
three point methods by some authors (jMaiid et al.l l2006l : 
llsmail et ailliooa '). 

Errors in a Semi-Lagrangian explicit-implicit integra- 
tion method include errors of the backward implicit integra- 
tion {0{At)), and errors from the interpolation: E{Ax) = 
0{Axk-i, Axk+i)- Therefore, the overall accuracy of the 
method is 



,n + l 



du 



— ^ ^ +o\ At 

At dt ' 



0(A<+^A<+^ 
At 



(46) 



An exact calculation of (47) can be found in 
l|Falcone and Ferrettilll998l '). Eq. 47 shows that the error is 
not monotonic with respect to At. Errors coming from a 
poor spatial interpolation could dominate. As At increases, 
the overall error decreases. If the first term 0( At) is insignif- 
icant, a further increase of n will not improve the overall 
accuracy. For the sake of simplicity, we consider negligible 
errors in the SPH spatial interpolation at the time step level 
n+l. This relevant argument is beyond the scope of this 
paper. On the other hand, if SPH spatial interpolation is in- 
adequate, even the explicit integrations do not improve the 
quality of results. Hence, if the solution is well-resolved in 
space with efficient spatial interpolations, the overall consis- 
tency (accuracy -I- stability + convergence) is solely due to 



the implicit integration method that, for a 3 FICS is well doc - 
umented in var i ous textbooks (for ex ample iFletcherl l|l99lh ; 
iLeVeaud (| 19921 . 12002D : iHirschI (jl993)). 



3.3 



Choice of a time step for an SPH 
Semi-Lagrangian explicit-implicit scheme 



A time step restriction is always necessary for time depen- 
dent calculations in computational fluid dynamics. Restric- 
tions are needed for mathematical stability in explicit calcu- 
lations. Instead they are necessary for accuracy considera- 
tions in implicit calculations. The Courant-Friedrichs-Lewy 
condition on the time step progression to solve PDE and 
ODE for explicit integration techniques offers a temporal ref- 
erence where numerical solutions are both stable and conver- 
gent with the mathematical solutions. Unfortunately, such 
a temporal referen ce is still debated for implicit integrati on 
of PDE and ODE l|Robertlll98l1 : IStaniforth fc C6telll99ll ). 
For SPH technique, the explicit time limiter is given by: 



At£ 



Cmmi=i,jv 



IV7 1-1 



1/2- 



, (47) 



which includes the Courant-Friedrichs-Lewy time lim- 
iter AtcFL. Vsig,ij is the signal transmission velocity be- 
tween close parti cles i and j within th e SPH spatial res- 
olutio n length h (|Monaghanlll985l . 1 1991 1 19971 : IWhitehurstI 
Il995l l. also including the sound velocity Cai, while \a\i is the 
full acceleration for the ith SPH particle. C is a number of 
the order of 0.2 - 0.5. 

Since in a correct Free Lagrangian particle fluid dynam- 
ics particles cannot interpenetrate with each other, another 
time step limit is the "kinetic" value: 



_h_ 
\f~ 



1/2 



(48) 



where | / 1 i | is the force per unit mass only due to gravi- 
tational and to pressure terms. However, such a longer lim- 
iter, could be correct for a Semi-Lagrangian implicit SPH 
hydrodynamics only for weak shock flows. To overcome such 
a difficulty, we consider the limiter: 

At, = (AtspHAtky/^ , (49) 

which corresponds to the geometric mean of the pre- 
vious two. This choice both largely prevents the numerical 
collapse of the time step to be used, and gives a good hydro- 
dynamics far from the risk of getting wrong solution for su- 
personic flows that occur when At ~ Atfe. Of course, for very 
extreme situations, even this time step limiter cannot be cor- 
rectly used, because of the absence of a temporal top limit 
for implicit integrations techniques. The ratio Atk / AtspH 
is always ^ 1. However, whenever Atk/ AtspH > 400 (e.g. 
At, > 20Atsp_ff), some consistent deviations from the cor- 
rect solutions are recorded. These deviations are much better 
confined in so far as Atk/ AtspH ^ 10. 

Hence, the choice of a time step for a Semi-Lagrangian 
explicit-implicit or for an implicit scheme is totally arbitrary, 
in principle, because it does not exist a theory predicting a 
time step limiter. For methods using implicit integration, it 



could also be possible to adopt t" 



t" 



At (X At SPH, 



whose parameter of proportionality should be found case 
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by case. Formulation (50) automatically simply correlates 
the implicit At to AtspH and to Atk- It is true that some- 
times we also could implicitly work using a larger At, but 
some cautions are needed, because accuracy of solution could 
surely be compromised. 

Given these definitions, the time step adopted for the 
Semi-Lagrangian explicit-implicit scheme discussed in §3.2 
is: -t"=At^ Ati - AtspH, being AtspH used in 
the first explicit step. Whenever At = 0, only the explicit 
technique is applied. 

Throughout the rest of the paper, we adopted C = 0.25 
in the evaluation of AtspH, Atk- A larger C value is always 
possible, up to C = 0.5 without relevant instabilities in the 
explicit solution. This implies that the longer AtspH will 
involve longer At;, according to the choice (50). 



4 CRITICAL TESTS 

What is previously discussed shows that the system of equa- 
tions of hydrodynamics in its Lagrangian form (eqs. 1-3) 
should be algebraically manipulated in a form similar to the 
finite difference methods, conserving the Lagrangian formal- 
ism. 

To verify whether the Semi-Lagrangian explicit-implicit 
approach is correct, discontinuity in the fiow represents crit- 
ical tests especially at the discontinuity front. The case of 
fiat translational fiows, or the propagation of linear sound 
waves, does not represent critical tests. In these cases, ei- 
ther trivial simplifications of the non-linear Euler system of 
equation, or negligible discrepancies in the fiow within SPH 
spatial resolution length h, detectable on scale length longer 
than h, do not produce sensible differences in the terms in- 
cluded in the continuity, energy and momentum equations 
(§3.2). To this purpose, we present both a SPH ID and a 
2D tests in the case of a blast wave and a SPH 2D test 
regarding a 2D shockless radial viscous transport in an an- 
nulus ring. The 2D modelling is also interesting to check the 
correctness and the consistency of the arbitrary statistical 
criterion of producing the integrated numerical solution by 
the average of as many integrals as many are the dimensions 
of the flow (§3.2). At the same time, these simulations allow 
us to verify whether the choice for the time step (§3.3) is 
free of criticisms. What is important is that SPH explicit 
and SPH Semi-Lagrangian explicit-implicit results of inte- 
gration strongly compare with each other. 

Even though a variable smoothing resolution length 
could also be used in SPH codes, the inadequacy of this 
approach in particular as for the radial transport whenever 
free boundary edge s characterize the fiow, r ecently discussed 
l|Lanzafamel l2010bl : iLanzafame et al.ll201ll ). suggests us in 
using a constant h. However, the adoption of an implicit in- 
tegration scheme to SPH is not correlated to its resolution 
length variability or constancy. 

Being the Semi-Lagrangian explicit-implicit integration 
involving a more laborious algorithm spending « 10% more 
time for each iteration, compared with the simpler explicit 
integration algorithm, there are not practical advantages 
in using the Semi-Lagrangian implicit integration mecha- 
nism in so far as Ati « 1.4 — l.^AtspH ■ However, even 
examples where the time step gain ratio is restricted only 
to 1 < Ati / AtspH ^ 2 are useful to understand the 



good and the bad of the implicit algorithm. Technic ally 
1 < Ati /At SPH ^ 6 l|Robertlll969l : [Robert et al.lll972l ) to 
avoid a relevant degradation in the accuracy of solution. 
This means that in so far 1 < Ati/ AtspH ^ 6 is respected, 
tests are meaningful. Otherwise some unexplored side of nu- 
merical techniques is faced. For this reason, tests here pro- 
posed are not so stressed out for working with a ratio of 
time steps beyond the limit of 6. A risk to be faced for 
computational research scopes is the adoption of a Semi- 
Lagrangian approach even if Ati/ AtspH ~ 5 — 10, betting 
that the accuracy in the solutions is quite goo d and that the 
larger numerical dissipation is limited (Robe7tl ll98lf ). In this 
case, the reduction of time in calculations is « 2 — 3 times 
compared with the explicit calculations. However, some- 
times, the gain-effectiveness of the Semi-Lagrangian tech- 
nique depends o n the context of the fluid dynamics problem 
(Bart ello and Th omas 1 996), as well a s on the spatial dis- 
cretization ( Falcone and Ferrettilll99^ ). 



4.1 ID and 2D blast wave tests 

The behaviour of shock fronts moving in the prevailing fiow 
is analytical ly described by the Rankine-Hugoniot "jump 
conditions" (lLeVeauelll99a |2002| : lHirschlll997j : lToroHl999l : 
lBatcheloill2000l '). These conditions are obtained by spatially 
integrating the hyperbolic Euler equations across the dis- 
continuity between the two fiow regimes left-right {I — r) in 
their Eulerian formalism. In ID, such equations are: 



d 



dx 
d_ 
dx 
d 



{pv) 

{pv^ + p) 
[pv{E + p/p)], 



(50) 
(51) 
(52) 



(53) 



dp 
dt 

dpv 
~dt 

dpE 
dt dx ^ 

where E — /2 + e, whose conservative analytical form 
can be synthesized as: 

dw d 
dt dx ■' 

Whenever in a shocktube the ratios pi/pr ~ ^i/^r ^ 1 
(and consequently p\/ pi = 1), and vi = Vr = 0, such a 
discontinuity is called a "blast wave". Fig. 2 shows a com- 
parison of both explicit and implicit SPH results together 
with the so called analytical solution. In the SPH blast wave 
test here considered, Pi/pi = ei/e2 ~ 10^, while the SPH 
particle resolution length is /i = 5 • 10~^. The whole compu- 
tational domain is built up with 2001 particles from X = 
to X = 100, whose mass is different, according to the shock 
initial position. At time T = all particles are motionless 
and the adiabatic index 7 = 5/3, while the ratios p\/ pi — 1. 
The first 5 and the last 5 particles of the ID computational 
domain, keep fixed positions and do not move. The choice 
of the final computational time is totally arbitrary, since the 
shock progresses in time. Initial values oi v, p and e at time 
T = are shown at the left-right edges of each plot of the 
same Fig. 2. 

Fig. 2 shows that both SPH results globally compare 
with each other and that they also compare with the ana- 
lytical solution. Both SPH (explicit and implicit) results are 
in a good comparison with the analytical solution. Discrep- 
ancies involve only 4 particles at most, with the exception 
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Figure 2. ID blast wave test regarding the comparison of both explicit (left size) and implicit (right size) SPH results to analytical 
(solid line) results. Density p, thermal energy e, pressure p and velocity v are plotted at time T = 5. Density and thermal energy of 
particles initially at rest at time T = refer to values plotted at the two edges for each plot. The initial velocity is zero throughout. 



of numerical solutions corresponding to analytical vertical 
profiles regarding thermal energy where, for both SPH solu- 
tions, discrepancies are larger. As Fig. 2 clearly shows, both 
the SPH numerical solutions suffer from some small well 
known i nstabilities, e specially in the proxim ity of disconti- 
nuities iGibbd 1 18981 . 1 18991 : iMonaghanI 1 19971 ) as far as the 
velocity profile is concerned. Such effect comes out when- 
ever a spatial high resolution is working together with an 
explicit handling of dissipation through an artificial viscos- 
ity damping to solve the Riemann problem of flow disconti- 
nuities. A low spatial resolution hides this effect because of 



the higher artificial damping due to a higher particle resolu- 
tion length h (eqs. 17-18). The higher the spatial resolution 
(the smaller h), the higher the "blimp" instabilities. More- 
over, in SPH, even the choice of the arbitrary parameters 
asPH and Psph should be linked to the specific physical 
problem. It is true that it is possible to slightly modify the 
artificial viscosity strength if only weak shocks appear in the 
problem at hand, thus, in general it is neither desirable nor 
necessary to tune the aspH and Psph parameters to a spe- 
cific problem. However, even asPH = 0.04 and Psph = 
were also used (|Meglicki et al.ll 19931 ) for turbulence develop- 
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Figure 3. Examples of density wave propagation for the 2D blast 
wave tests. The density wave front develops from X = at time 
T = 0. Only plots at times T = 1,3,5 are represented. Density 
isocontour maps in 64 greytones are shown in each side. 



Figure 4. Ati/ AtspH and relative errors in the total energy 
AE/E per unit mass, from the total energy E conservation, re- 
ported as a function of time T for the 2D blast wave for both 
explicit and explicit-implicit integrations. Notice that for the sake 
of simplicity of representation, these errors are multiplied by 10"^.. 



ment. Another way f or reduction of instabili ties is possible 
ijLanzafamd I2010al l3: iLanzafame et all l201lh if the damp- 
ing is strictly locally physical, using eqs. (23, 24) instead 
of eq. (19) for the perfect gases EoS. An effective reduc- 
tion of "blimp" effects is obtained, especially for ID blast 
waves, where strong discontinuities in the flow deeply affect 
the SPH numerical solution producing intrinsic numerical 
instabilities close to the propagating discontinuities. 

To further check the consistency and the validity of the 
Semi-Lagrangian explicit-implicit integration method, we 
also performed both explicit and explicit-implicit 2D blast 
wave tests on the basis of the same initial physical param- 
eters, considering an initial ordered disposition of particles, 
adopting the same particle resolution length h and the same 
Ax = h = 0.05 and Ay — h — 0.05 particle separation. The 
four edges of the box at Xmi„ = —20, x^ax = 20, ymi„ = 0, 
Vmax = 20 are formed of three equally spaced lines of always 
motionless particles. As in the previous example, the blast 
wave starts from a; = 0. The temporal progress of the den- 
sity wave is shown in Fig. 3 for both the integration tech- 
niques. For the sake of simplicity, only grey tone contours 
are shown at time T — 1,3 and 5. The agreement looks 
like quite good. However, as we discussed in §3.1, despite 
stable, the implicit integration techniques involves greater 
errors compared with those made in explicit integrations. 
Fig. 4 shows the relative error for both integration schemes 
on the total energy AE/E, where it is clearly shown that 
in the Semi-Lagrangian explicit-implicit integration scheme, 
relative errors on the total energy per unit mass per parti- 
cle E — w^/2 -f e are larger, being of the order of ~ 2 • 10"'' 
against ~ 2 -10"^, relative to the explicit integration scheme, 
as can be evaluated by the two pictures shown in this figure. 
However, what is relevant is that in the Semi-Lagrangian 



explicit-implicit integration method such error does not in- 
crease in time compared with that relative to the explicit 
integration technique. The most of the relative error is made 
at the first stages of integrations in both cases, after that it 
looks like to be stationary. 

Apart from the accordance between analytical and nu- 
merical (explicit and implicit) results, what is remarkable 
is the fact that implicit results are obtained in about half 
of cpu time in implicit SPH approach, compared with the 
traditional explicit SPH, adopting the time step criterion 
expressed by eq. (46). The final time configuration at time 
T = 5 is an evolved configuration well beyond the initial 
instants from T = 0. At the final time configuration at time 
T = 5 the time step gain ratio Ati/ AtspH ^ 1.6. Since 
the beginning of a blast wave simulation, a time step gain 
ratio ~ 2 is quickly achieved as an order of magnitude. At 
the beginning of simulations, the time stepping is principally 
governed by the thermal contribution thanks to the sudden 
increase in the velocity. Time by time, the decrease of v and 
the relative decrease in At^, set a larger and larger Ati (de- 
rived by Atk) with respect to AtspH- This implies that in 
the particular case of a detonation, the criterion adopted on 
the time step progress {Ati in §3.3) shows accordance with 
analytical results without any difficulty, being Ati/ AtspH, 
well within the factor of 6 discussed in §1. Notice that this 
accordance is not verified at the 1st time step, because an 
explicit integration scheme must be applied. 

4.2 2D shockless isothermal radial viscous 
transport in an annulus ring 

Theory on 2D shockless radial transport in a Keplerian an - 
nulus ring in a gravitational potential well (jPringlel [l98lh 



Implicit integrations in SPH modelling 13 




r r 



Figure 5. Comparison of radial density distributions of a viscous 
annulus ring for both explicit (left size) and implicit (right size) 
SPH results in the Shakura and Sunyaev formulation using agg = 
0.4. Viscous time 9 for each distribution is also reported. 

predicts that, from the Green function, the solution of the 
initial Keplerian mass distribution at time t = for E is: 

E(r-,t = 0) =m(5(r- To )/27rro (54) 

for a ring of mass m, with an initial radius ro, The 
solution, at time t, in terms of dimensionless radius x = r/ro 
and viscous time 9 — 12vtr^^ is 

T.(x,t) = (m/7rro~^)6»~'x"'''^exp[-(lW)/6l]/i/4(2x/6»),(55) 

where I^/^ is the modified Bessel function. The action 
of viscosity is to spread out the entire annulus ring toward 
a disc structure transporting most of the low angular mo- 
mentum mass toward the centre of the potential well and 
transporting a smaller fraction of high angular momentum 
mass toward the empty external space. 

For practical computational purposes, since it is impos- 
sible to reproduce a delta Dirac function at time t = for 
the initial mass distribution, it is necessary to start numer- 
ical calculations from an initial mass distribution relative 
to a small 9 value. In our exampl e, we choose 9 = 0.01 7, a 
value comparable to that used bv lSpeith fc Kiev! (|2003l ) for 
a radial viscous transport similar test. 

The initial setting is that of a 2D annulus ring includ- 
ing 40000 particles whose h = 0.09 forming a ring whose 
radius ro = 1 around a star having mass M = 1. The initial 
radial spread &t 9 = 0.017 is Ar = 1. Throughout the simu- 
lations de/dt = 0, that means that the entire configurations 
is permanently isothermal. 7 = 5/3 and the unifirm sound 
velocity Cs = 5 to have low Mach numbers from the begin- 
ning of the simul ations. In realit y , a va lue of Cs — 0.05, like 
that adopted by ISpeith fc kIctI (|2003D would not address 
any time stepping advantage to the implicit integrations be- 
cause Ati ~ At SPH in isothermal conditions, preventing 
any flow heating. For a dimensional arbitrary initial setting, 
we consider a central star having Mo — IMq, while the 
dimensional annulus ring radius Ro = 10^^ cm. These are 
also normalization factors for masses and lengths for a di- 
mensionless handling of equations, as previously discussed 
in §3.2. Vo = 2-k(GMo/Ro)^^^ is the normalization factor for 
velocities, so that the orbital period is normalized to 1. 

Figure 5 compares results on the radial distribution 
of the 2D mass density E for the viscous annulus ring at 
various viscous times 9 both for SPH explicit integrations 
a nd for Semi-La g rangian explicit-impl i cit in tegrations for 
a lShakural l|l972l ): IShakura fc Sunvae^ (|l973l ) viscosity pa- 
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Figure 6. Comparison of relative errors for the total energy E 
per unit mass and for the specific angular momentum r^uj, for 
the radial viscous transport in an annulus ring as a function of 
the viscous time 8. Explicit and Semi-Lagrangian explicit-implicit 
relative errors are compared. 








Figure 7. Implicit to explicit time step Ai; / Atspu reported as 
a function of the viscous time 6 for the viscous spreading of the 
annulus ring. 

rameter ass = 0.4. In figure 5, these two pictures clearly 
show that, T,{r,9) for implicit integration suffers of the cu- 
mulation of both explicit and iterated implicit dissipation 
viscosity due to the effect of a further intrinsic viscosity al- 
gebraically developed by repeated implicit iterations. Even 
though fortunately limited, relative errors in figure 6 for the 
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implicit integration scheme on the total energy per unit mass 
E = /2 + e — l/r and on the specific angular momentum, 
as a function of viscous time 0, are moderately larger com- 
pared with those relative to the explicit integration scheme. 
This is not a surprise, since it is previously discussed why 
errors in the implicit integration techniques are larger. The 
moderately larger inaccuracy is however compensated by a 
larger stability in results and by a shorter computational 
time. In the example, here discussed, we yielded implicit re- 
sults after a total time ~ 1/2 shorter than the total time to 
conclude explicit calculations. In Fig. 7, it is shown the time 
step gain ratio Ati/AtspH as a function of the viscous time 
8. At the beginnning of the simulation the explicit AtspH 
is substantially due to the thermal contribution, so that the 
ratio between the two time steps is at its maximum, being 
the tangential kinematics substantially Keplerian and being 
the annulus ring still at its initial configuration. Once time 
progresses, At^ becomes shorter as a consequence of the in- 
ner edge Keplerian kinematics, so that differences between 
the explicit and the implicit time stepping reduce. In this 
example, the decrease of the kinematic At^ time by time 
caused by the inner edge contraction toward shorter r, re- 
duces differences between the thermal explicit AtspH to the 
implicit Ati. 

Results of this test clearly show the good and the 
bad side of a Semi-Lagrangian explicit-implicit integration 
techique compared with a simpler explicit integration tech- 
nique in the worst situation when in the case of isother- 
mal flows, there is not any advantage in using an implicit 
time stepping. Solution are stable and accurate, in so far 
as Ati/ AtspH is not too large as previously discussed, but 
results are affected by a larger numerical dissipation devel- 
oped by the iterative implicit integration procedure. It is 
this larger dissipation that ensures a better stability of so- 
lutions at the price of a moderate inaccuracy of solution. 
Despite this example shows how the implicit integration is 
not competitive to the explicit technique in these situations, 
results here shown illustrate how much is the difference in 
explicit and implicit results in the worst competitivity con- 
ditions to the implicit integration scheme. The presence of 
coUisional shocks, here prevented, would decrease AtspH, 
better stressing the difference between AtspH and Ati. 



5 SIMULATIONS OF A 3D ACCRETION DISC 
AROUND A MBH IN A CLOSE BINARY 

We compare results of both inviscid and viscous station- 
ary disc structures performing SPH simulations whose inte- 
gration schemes are both explicit and implicit in low com- 
pressibility (7 — 1.3) with the aim of getting a physically 
well-bound accretion disc around a MBH a close binary. 
Previous preliminary results on this them e were published 
(|Lanzafame fc Belvedere! l200ll . |2002| . |2005| ) both in 2D and 
in 3D. 

The characteristics of the binary system are determined 
by the masses of the MBH and of its companion star and 
their separation. We chose to model a system in which the 
mass Ml of the primary MBH and the mass M2 of the sec- 
ondary subgiant star are equal to 32M0 and IM©, respec- 
tively and their mutual separation is di2 = 10* Km. The 
primary's potential well is totally empty at the beginning 
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Figure 8. XY plots of the injected particle stream for the non 
viscous (ogs = 0) and the two disc models {ctss = 0.1 and 
agg = 0.4) at time t = I. The central MBH position is also 
shown. Notice that the Y positions are shifted for agg = and 
for agg = 0.4. 



of each simulation at time T = 0. The injection gas veloc- 
ity at LI is fixed to Vinj — 130 Km s"^ while the mjec- 
tion gas temperature at LI is fixed to To — 10* K, taking 
into account, as a first approximation, the radiative heat- 
ing of the secondary surface due to radiation coming from 
the disc. Gas compressibility is fixed by the adiabatic index 
7 = 1.3. Super s onic kinema t ic conditions at LI are di scussed 
in iLanzafamd (|2008l . l2009l l: iLanzafame et al.l l|2006l ), espe- 
cially when active phases of CB's are considered. However, 
results of this paper are to be considered as a useful test to 
check whether disc structures (viscous and non) show the 
expected behaviour. The reference frame is centred on the 
primary compact star, whose rotational period, normalized 
to 2n, coincides with the orbital period of the binary sys- 
tem. This explains why in the momentum equation (eq. 2), 
we also include the Coriolis and the centrifugal accelerations. 

Pressure, density, temperature and velocity are six un- 
knowns to be found. Therefore we solve the continuity, mo- 
mentum, energy, and EoS equations. In order to make our 
equations dimensionless, we adopt the following normaliza- 
tion factors: M — Mi + M2 for masses, di2 — W^^ cm 
for lengths, Vo = (G(Mi -I- Af2)/di2)^''^ for speeds, so that 
the orbital period is normalized to 2tv, po — g cm~^ 

for the density, po = povl dyn cm~^ for pressure, vl for 
thermal energy per unit mass and To = (7 — 1)110 Wp Kg^ 
for temperature, where rrip is the proton mass and Kp is 
the Boltzman constant. The adopted Kernel smoothing res- 
olution length is ft = 5 ■ 10"'^ throughout. The geometric 
domain, including disc particles, is a sphere of radius 1, cen- 
tred on the primary MBH. The rotating reference frame is 
centred on the compact primary and its rotational period 
equals the orbital one. We simulated the physical conditions 
at the inner and at the outer edges as follows: 
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Figure 9. XZ plots of the injected particle stream as in Fig. 8. 

a) inner edge: 

the free inflow condition is realized by eliminating parti- 
cles flowing inside the sphere of radius 10~ , centred on the 
MBH. Although disc structure and dynamics are altered 
near the inner edge, these alterations are relatively small 
because they are balanced by a high particle concentration 
close to the inner edge in supersonic injection models. 

b) outer edge: 

the injection of "new" particles from LI toward the in- 
terior of the primary Roche Lobe is simulated by generat- 
ing them in fixed points, called "injectors", symmetrically 
placed within an angle having LI as a vertex and an aper- 
ture of ~ 57°. Normally, as adopte d since our fir st paper 
on SPH accretion disc in CB (.Molteni et al.lll991^ . the ra- 
dial elongation of the whole group of injectors is ~ IQh. 
The initial injection particle velocity is radial with respect 
to LI. In order to simulate a constant and smooth gas injec- 
tion, a "new" particle is generated in the injectors when- 
ever "old" particles leave an injector free, inside a small 
sphere with radius h, centred on the injector itself. Parti- 
cle masses are determined by the assumed local density at 
the inner Lagrangian point LI: pLi = 10~®gi cm~^ (as typ- 
ical stellar atmospheric value for the secondary star), equal 
to m = pLi{hdi2f /{Ml + Ma). 

The adoption of super sonic mass transfer co nditions 
from LI is fully discussed in iLanzafamd ||2008| . |2009D . where 
disc instabilities, responsible for disc active phases of CB are 
discussed in the light of local thermodynamics. Whenever a 
relevant discrepancy exists in the mass density across the 
inner Lagrangian point LI between the two stellar Roche 
lobes, a supersonic mass transfer occurs as a consequence of 
the moment um fiux conservation . The same result can also 
be obtained (iLubow fc iSh3[l97i) by considering either the 
restricted problem of three bodies in terms of the Jacobi con- 
stant or the Bernoulli's theorem. Moreover, and this is the 
most important thing, we need to compare 3D disc mod- 



els where AtspH and Atk could be significantly difi^erent. 
This condition is searched for since the injection conditions 
at LI, favouring violent collisions among low compressibility 
(7 — 1.3) particles moving around a MBH, rising up particle 
heating at expense of the kinetic -I- gravitational energies. 
These conditions make much smaller the thermal contribu- 
tion in At SPH evaluation compared with the kinetic con- 
tribution throughout the disc, especially when viscous heat- 
ing or other forms of heating are considered. Thus when 
Ati/ AtspH ^> 1, a significant deviation in the implicit so- 
lutions would afi^ect the whole result, up to compromising 
the numerical stability. A sensible reduction in AtspH even 
happens whenever other forms of signal transmission veloc- 
ity (e.g. the Alfen speed) are also taken into account. This 
makes the inner disc region, at its inner edge, quite criti- 
cal for a low compressibility highly viscous accretion disc 
around a MBH, because all these characteristics contribute 
to a sharp decrease of the explicit AtspH- In particular a 
supersonic mass viscous fiow transfer at Li and a central 
MBH determine an initial mechanical energy large eonugh 
to be converted into heat with high 7 to yield very low Mach 
numbers at the disc inner edge. 

The first viscosity coefficient (eq. 10) in the viscous disc 
mode ls is related to t he Shakura and Sunyaev p arametriza- 
tion (|Shakural [l973 : IShakura fc SunvaevI [l97a ) as: rjv — 
passCsH , being H the local disc thickness. The second vis- 
cosity coefficient, related to the bulk viscosity is not taken 
into account for the sake of simplicity. Two viscous disc mod- 
els are considered, whose ass = 0.1 and whose ass ~ 0.4. 
These values are in accordance with the typical and with 
the maxim um ass com patible with both astrophysical ob- 
servati ons (iKing et al.l l2"007i'l and w ith laboratory experi- 
ments l|Abolmasov fc Shakurall2009l ) . These disc models are 
also compared with the non viscous disc model counterpart, 
taken as a reference model. 

Figg. 8 and 9 show three XY and XZ plots of the in- 
jected stream from LI at time t = 1 for the three disc mod- 
els. The greater compactness of the injected flow of parti- 
cles is visible with the increase of the viscosity parameter 
as a result of the sticking viscous effect. The large stream 
geometric spread is mainly due to the low compressibility 
here adopted (7 = 1.3), while the wide initial circulariza- 
tion radius is explained by the high angular momentum at 
LI of the injected flow. These initial kinematic conditions 
affect the whole disc structure and kinematics throughout 
the simulations. Indeed, the injected gas stream yields an 
impenetrable boundary of the outer disc edge itself and a 
significant fraction of the disc's ejection flow comes from 
this side of the outer disc edge. Besides, a wide spray of di- 
rect cold injected flow above and below the mean disc plane 
also plays a role. 

Figs. 10, 11 and 12 show a comparison of explicit and 
implicit SPH disc structures for the non viscous and the 
viscous models whose ass =0.1 and 0.4, respectively. All 
explicit and implicit structures impressively compare with 
each other, as well as their injection, ejection and accre- 
tion rates, whose values are in the order of ~ 10^^ g , 
6.5 ■ 10^° g s"^ and 3.5 • 10^° g respectively for the 

inviscid disc models, of ~ 7.8 ■ 10^° g s'S 3.6 ■ 10^° g s"^ 
and 4.2 ■ 10^° g s~^, respectively for the ass ~ 0.1 vis- 
cid disc models, of ~ 5.9 ■ 10^° g s'S 2.6 ■ 10^° g s"^ and 
3.3 ■ 10^° g , respectively for the ass ~ 0.4 viscid disc 
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Figure 10. XY plots of 64 grcytoncs density p isocontours of the 
non viscous 3D SPH disc modelling in a microquasar. Explicit 
SPH results compare with SPH implicit ones. 



Figure 11. XY plots of 64 grcytoncs density p isocontours of the 
oiSS = 0.1 viscous 3D disc modelling in a microquasar. Explicit 
SPH results compare with SPH implicit ones. 



models. Some limited differences appear in the comparison 
between the implicit to the explicit disc structures in Figg. 
10 and 11, due to the fact that the implicit structures suf- 
fer of a relative additional intrinsic viscosity, because of the 
cumulation of numerical dissipations in the implicit integra- 
tion iterative procedure. These differences are reduced and 
become negligible in the physically more viscous disc struc- 
ture of Fig. 12, where both the explicit and the implicit disc 
structures better compare with each other. In this case the 
higher physical viscosity makes secondary the effect of other 
artificial and intrinsic numerical dissipations. 

The total number of disc particles are of the order of 
137613, 143013 and 125203 particles, respectively in steady 
state conditions when the mass of the disc is statistically un- 
changed for explicit calculations after 522100, 5246800 and 
7153700 time steps for ass = 0, 0.1 and 0.4, respectively. 
Instead, for implicit calculations, we count 135112, 139717 
and 123150 particles, respectively after 148300, 178250 and 
248150 time steps, respectively for the same ass- As it is 
evident, the total number of particles compare with each 



other for the same ass characterizing the disc model, but 
the total number of time steps for implicit integrations is 
clearly much less than that relative to explicit integrations. 

The durations of each integration process for the im- 
plicit scheme is of course longer taking into account of the 
iterative cycles. The total number of disc particles for im- 
plicit integration schemes are about 1 — 2% systematically 
lower than those relative to explicit integration schemes. Of 
course these are very little difference among explicit-implicit 
results. This result is explained by the fact that, being the 
mass transfer process from LI a discrete process simulated 
by a particle generation, the longer implicit time step in- 
volves a mass transfer rate about 1 — 2% slower in the im- 
plicit cases. Consequently, also the accretion and the ejec- 
tion rates for implicit integration disc modelling are propor- 
tionally lower than those relative to the explicit integration 
disc models. Being these differences so tiny, they are not de- 
tectable by the direct counting of particles on short or on 
medium size time intervals, refereed to the orbital period. 
However, this very little difference in the injection, ejec- 
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Figure 12. XY plots of 64 greytones density p isocontours of the 
Q5S = 0.4 viscous 3D disc modelling in a microquasar. Explicit 
SPH results compare with SPH implicit ones. 



tion, accretion rates, does not involve consistent difTerences 
among the explicit - implicit disc structures and dynamics, 
whose azimuthal and radial profile of density look like very 
similar, a part some graphical contrast effects in the grey 
tones. 

Even though disc models refer to a low compressibil- 
ity regime (7 = 1.3), the primary's MBH Roche lobe is 
large and deep enough to favour a disc consistency even 
in the non viscous disc model. This is a first result, different 
from those relative to low mass binaries, where SPH mod- 
els yield sc arcel y po pulated low compressibility non vi scous 
structures (|Molteni et al.lll99ll : iLanzafame et al.lll992l ). Al- 
though finalized to a different strategy, a statistically signif- 
ica nt 2D structure of the disc arou nd a MBH were obtained 
by fcanzafa me fc BelvederlbOOSl ). 

The colhsional push exerted by the flow coming from LI 
on the outer edge of the disc yields an effective perturbation 
generating a global disc's elliptic geometry in the viscous 
cases. As a consequence, spiral density patterns character- 
ize SPH viscous disc models around MBH, as well depicted 
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Figure 13. XY plots of the subsonic spiral kinematics ot the 
agg = 0.1 viscous 3D disc in a microquasar. Selection in the 
radial Mach number are shown in each panel. 

in Figg. 11 and 12. Instead, possible spiral patterns are not 
well resolved, or they do not well develop in the SPH non 
viscous model because of the inadequacy of dissipation able 
to numerically resolve shock fronts in the low compressibil- 
ity flow and because the stream flow from LI is too sparse 
in order to exert an effective localized push at the disc's 
outer edge. In these two last figures in particular, relative 
to the two physically viscous models, two main spiral pat- 
terns in the density are evident, even if the glimmer of the 
appearance of a third one is also shown. These facts are 
recorded in 3D modelling. On the contrary, spiral structures 
around MBH, as well as spiral shocks in the disc's radial 
flow, are normall y seen in the 2D non viscous lo w compress- 
ibility modelling (|Lanzafame fc Belvederell2005l ) . In the case 
of low mass close binaries, SPH low compressibility viscous 
disc models are not able to show any spiral pattern because 
of several shortcomings in the method, especially when free 
edge conditions are considered l|Lanzafamell2010bl ). 

The better compactness of the entire elliptical viscous 
disc structure, not only yields disc's spiral density proflles, 
but contextually also a spiral radial flow kinematics. In these 
low compressibility viscous models, this is a subsonic kine- 
matics, as shown in Figg. 13 and 14 where a selection in the 
radial Mach number is made. The tendency should be to pro- 
duce spiral shocks which should appear either increasing the 
flow compressibility modelling (decreasing 7) (Chakrabart^ 
Il992l l. or altering the stellar mass ratio or the initial injec- 
tion kinematics in order to simulate a mass transfer flow 
from LI where the radial kinematics is more enhanced than 
the tangential one. The higher are the kinetic energy and the 
initial angular momentum, t he better is the coming out o f 
spiral structures and shocks (|Lanzafame et ai1l2000l . [2OOII ). 
This is a consolidated result that has emerged for low mass 
binary systems. However, in the case of a microquasar, where 
the primary component is a MBH, a large initial angular 
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Figure 14. XY plots of the subsonic spiral kinematics ot the 
c«SS = 0-4 viscous 3D disc in a microquasar. Selection in the 
radial Mach number Mj. are shown in each panel. 

momentum and large Coriolis and centrifugal terms favour 
too much an initial tangential kinematics at the cost of the 
radial one. 

A further comprehension of disc structure and kinemat- 
ics is shown in Figg. 15 and 16, where the radial distribution 
of the specific angular momentum and the radial distribu- 
tion of temperature are respectively shown on a logarithmic 
scale for the explicit integration results only because those 
relative to the implicit integration process are substantially 
the same. Two details come out from these pictures. The 
first one is that the inner disc regions are less populated of 
particles in the two viscous models. This effect, due to the 
enhanced radial viscous transport, makes the radial profile 
of these two distributions more "flat" compared with that 
relative to the non viscous distributions. The second partic- 
ular that appear from these two pictures is that the slopes 
of the two radial distributions are lesser than r^Q oc r^^^ 
and T oc r~^^* relative to the "standard disc model". The 
temperature radial profile is even flat for a large portion of 
the external part of the discs. This result is a consequence of 
the wide geometric opening of the injected flow coming from 
LI (Figg. 8-9), whore colder and higher radial flow in the 
disc bulk mix with hotter flows transported from the disc 
outer edge toward the central accretor. As a consequence, 
looking from the outer disc edge toward the centre of the 
disc, both the temperature radial increase and the specific 
angular momentum decrease are lower than T oc r"^^** and 
r^Q oc r^^^ laws, even though the viscous heating effect in 
the disc bulk and in its inner regions is remarkable especially 
for the ass = 0.4. 

From the pure numerical point of view, the so tight 
correspondence among explicit and implicit disc structures 
means that the radial viscous transport mechanisms are 
comparable despite the higher accumulation of dissipation 
and error propagation during the implicit integration it- 



o.tr~i — ' — I — I — I — I — I — I — I — I — ' — ' — 3 — ' — I — T — I — J — J — z 



-o.z r 




J;, ■ 'ji — J ] 1 — 1 \ 1 1 — a — J — J — ■! — J — ;i — a — ] — a — a — a — a — 



Figure 15. Radial distribution of the specific angular momentum 
in a logarithmic scale. 
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Figure 16. Radial distribution of temperature in a logarithmic 
scale. 



erative cycles. This is explained by the fact that the in- 
tegrated effect of dissipation, distributed on a At; up to 
10 — 20 times longer than Ats pu accumulated on 3 — 4 cy- 
cles corresponds to that relative to explicit integrations dis- 
tributed on the same time interval (At; = mAtspH, where 
m = Ati/AtsPH)- This non trivial numerical conclusion 
confirms that, at least for SPH, the chosen criterion for the 
number of iterative cycles works well. 
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6 CONCLUDING REMARKS 

From the numerical point of view, a successful Semi- 
Lagrangian explicit-implicit integration numerical scheme 
is here applied f or the SPH meth od that is a Free La- 
grangian scheme (|Whitehurstll 19951 ) . Comparison to results 
obtained working with explicit numerical schemes shows an 
impressive convergence of results. Traditionally, numerical 
schemes in such explicit-implicit approaches are stable and 
show a convergence when the adopted implici t time step is 
larger than the exp licit one by up to 6 times (|Robertll 19691 : 
[Robert et al.|[l973 ). In this study, correlating our implicit 
time Ati stepping to the explicit time stepping AtspH and 
to the kinematic time stepping Atk, we report consistent and 
stable results both for some ID and 2D critical tests involv- 
ing shocks and shockless physically radial viscous transport, 
as well as for 3D low compressibility accretion disc models 
around a MBH in a microquasar. This choice, in particular 
for the viscous modelling is motivated to stress the validity of 
the implicit technique we described. Consistency of results is 
recorded up to Ati/AtsPH = 15 as shown in Fig. 17 for the 
cess = 0.4 physically viscous accretion disc modelling, where 
Ati/ AtspH is reported as a function of time for the implicit 
results for the three simulated discs. This has the evident ad- 
vantage that working in such an explicit-implicit scheme we 
can obtain meaningful results in a much shorter time than 
working adopting the SPH explicit time stepping. In particu- 
lar we yield implicit results for the ass = 0.1 and ass = 0.4 
discs in two weeks of cpu time using a serial code on a sys- 
tem based on an AMD Opteron cpu instead of times of the 
order of 5 months on the same hardware working in explicit 
integration technique. Compared with other implicit tech- 
niques, this technique, working with a larger time computa- 
tional step than the AtspH, is competitive in so far as the 
integrated numerical value is obtained in a few iterations, 
otherwise the time involved in a higher number of iterations 
makes the techniques not attractive. The constraint for our 
Semi-Lagrangian approach we propose is that the number 
of implicit iterations should be high enough to have stable 
solutions, but low enough to have accuracy with a moderate 
dissipation and not too time consuming. However, a fast con- 
vergence of the integrated numerical solution belongs to the 
SPH philosophy because of the spatial smooth distribution 
of any physical property around each moving particle, where 
V j^A{r')W{r,r\h)dr' ~ A{r')WW{r,r' ,h)dr' . At 
the same time, the necessity to compute some further spa- 
tial gradients according to eqs. (45, 46) does not involve any 
further cpu time because these spatial gradients could be 
calculated contextually with other SPH quantities during 
the explicit step described in §3.2. 

On §3.2 we put an unanswered question on which im- 
plicit relaxed solution is obtained after numerous iterations 
especially when 10 < Ati/ AtspH < 20. This argument 
of pure numerical applied mathematics is worth to be the 
theme of a dedicated book on the stability and conver- 
gence of Semi-Lagrangian implicit numerical solutions of 
hyperbolic - parabolic systems of equations in fluid dy- 
namics and is well beyond the scope of this paper, whose 
message is that for SPH it is possible to calculate consis- 
tent Semi-Lagrangian implicit solutions even in the range 
10 < Ati/ AtspH < 20 working with 3 — 4 iterative cycles 
(the first one explicit). Hence, these results further shift the 




Figure 17. Ati/AtspH as a function of time t for the three 
simulated 3D accretion discs. The panel at the bottom refers to 
the ass = i^oii viscous model, agg = 0.1 (middle panel) as well 
as ass = 0-4 (top panel) results are also reported. 



limit of 6 claimed by 'Robert' (1969'); 'Robert et al.' (197^ 
and stiU debated Robert (1981); Staniforth & Cote (1991). 

From the astrophysical point of view, a well-bound ac- 
cretion disc appears in a microquasar even for a non viscous 
- low compressibility modelling. This is allowed because of 
the strong gravitational field of the central MBH prevail- 
ing on the repulsive pressure forces. In these low compress- 
ibility models the kinematic of the radial flow stays sub- 
sonic throughout even if the physical viscous contribution 
is also taken into account. For the physically viscous mod- 
els, two main spiral patterns in the density flow clearly ap- 
pear mainly as a consequence of the elliptical morphology of 
the entire disc structure and the disc's outer edge is better 
shaped. A possible third spiral profile can also appear. The 
radial disc structures all deviate from that of the disc stan- 
dard model because of the vertical inflow of cold gas coming 
from LI. 

The full consistency of the Semi-Lagrangian explicit- 
implicit integration technique applied to free Lagrangian 
methods is checked not only because of the technical tests 
here presented, but above all, because of the fully con- 
vergence of results between 3D SPH explicit and explicit- 
implicit accretion disc modelling, where a complex of phe- 
nomenology, both viscous and inviscid, are here shown, es- 
pecially when a strong gravitational field of a MBH deeply 
affect the local viscous thermodynamics producing a consis- 
tent difference between Ati and AtspH- This beyond the 
correspondence between our disc modelling and a real exist- 
ing astrophysical counterpart. Around a MBH, some physics 
of viscous disc models is stressed so much to determine, es- 
pecially at its inner edge those conditions decreasing the 
explicit AtspH to numerical values so small to stress out if 
a Ati/AtsPH > 6 could yield to SPH consistent solutions 
in a Semi-Lagrangian explicit-implicit integration. 
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A gravitational field, even if dominant against pres- 
sure forces is necessary to check if for large Atk/^tsPH the 
method correctly works. Moreover, in a 3D problem, rele- 
vant gravitational forces could be dominant against pres- 
sure, viscosity ones, etc., only along the radial direction. 
The opposite happens along the tangential and vertical di- 
rections. This as for the momentum equation only. In any 
case, the consistency of solutions should also take into ac- 
count of the continuity and of the thermal energy equations 
not including any external field. Since the complete solu- 
tion must also take into account of these contributions, we 
deduce that the comparison of explicit and explicit-implicit 
solutions is a success of the new scheme. 
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